Week 3 Lab - Maps

#########################################################################################
################################# GEOGRAPHIC MAPS IN R ##################################
#########################################################################################

## OUTLINE ##

# 1.  National and Subnational Maps Using Tigris
# 2.  World Maps Using Rnaturalearth
# 3.  National and Subnational Maps Using Tigris and Tidycensus for Analyzing Census Data

# Install these packages (if you haven't already):
# install.packages("tigris")
# install.packages("tidycensus")
# install.packages("rnaturalearth")
# install.packages("ggthemes")
# install.packages("tidycensus")
# install.packages("ggrepel")
# install.packages("socviz")
# install.packages("mapview")

library(tidyverse)
library(tigris)
library(rnaturalearth)
library(tidycensus)
library(ggthemes)
library(socviz)
library(poliscidata)
library(scales)
library(haven)
library(ggrepel)
library(mapview)

#########################################################################################
################### 1.  National and Subnational Maps Using Tigris ######################
#########################################################################################

# Go-to source for tigris: Analyzing US Census Data: Methods, Maps, and Models in R, 
# by Kyle Walker. https://walker-data.com/census-r/index.html
# See in particular Ch. 5 for basic information on tigris
# This book is the same go-to source in Part 2 on "tidycensus"

# Q:   What is a CHOROPLETH MAP? 

# Start simple: Use tigris to get national map with states (take out PR)
st <- states(cb = TRUE, resolution = "20m", progress_bar=FALSE) |>
  filter(NAME != "Puerto Rico") 

# View

# Check "class" of this data/maps object; you'll see it's "special features" (sf) and 
# a data frame sf = contains geographical information for mapping.
class(st)

## [1] "sf"         "data.frame"

# You can sort by state name if you want. This is not necessary, however.
st <- arrange(st, NAME)

# Quick map (notice how much simpler the syntax is using "sf" method)
ggplot(st) + geom_sf()

# What happened?! :-) 

# Add "shift_geometry" at the end; gives better look, and puts AK and HI on bottom
st <- states(cb = TRUE, resolution = "20m", progress_bar=FALSE) |>
  filter(NAME != "Puerto Rico") |>
  shift_geometry()

# Quick map: Geometry here is "sf", which means ggplot recognizes this is a 
# "special features" geographic maps object.
ggplot(st) + 
  geom_sf() 

# Quick map - theme_void to remove gray background
ggplot(st) + 
  geom_sf() + 
  theme_void()

# Change colors: Check out the color palette file in the Week 4 data session folder
ggplot(st) + 
  geom_sf(fill="dodgerblue") + 
  theme_void()

# Change state border line colors
ggplot(st) + 
  geom_sf(fill="dodgerblue", color="white") + 
  theme_void()

# Change state border line width - thick
ggplot(st) + 
  geom_sf(fill="dodgerblue", color="white", linewidth=.8) + 
  theme_void()

# Change state border line width - thin
ggplot(st) + 
  geom_sf(fill="dodgerblue", color="white", linewidth=.1) + 
  theme_void()

# Let's start with a very simple choropleth map - 2020 election, Trump v. Biden
# Import 2020 state election data
v2020 <- read_csv("us_vote_2020.csv")

# View

# Our task in generating a choropleth map: 
#   We need to merge our election data (v2020) with our "sf" maps object, "st"
#   See Ch. 3 in ModernDive, section 3.7 (using the "join" function) 

#   First, match state names across datasets; use the "mutate" function, which generates
#   a new variable in the data object. This sets up for the merge, which requires you to
#   have **one** identification variable in common between the two datasets. In this case,
#   that's going to be "NAME."
v2020 <- v2020 |> mutate(NAME=state)

# Note how I'm using piping; the above command would be identical to: 
#   v2020 <- mutate(v2020, NAME=state)
#   Remember, when you use piping, you call on the dataset before the pipe operator
#   and then in subsequent commands, you don't need to call in the data.

# Merge our v2020 data object with our st object containing geographic information using
# "left_join." Notice the syntax there. Put the "master" data object first (i.e., st) and 
# the data object that you're merging into the master object second.
# "left_join" automatically merges by the variable that is common between the two
# datasets. Again, in this example, that's "NAME."
stv20 <- left_join(st, v2020)

# Map using "theme_void"
ggplot(stv20, aes(fill = called)) + 
  geom_sf(color = "gray90", linewidth=.2) + 
  scale_fill_manual(values = c("blue", "red"), labels=c("Biden", "Trump")) + 
  theme_void() + 
  coord_sf(expand=FALSE) +
  labs(fill = "",
       title = " 2020 US presidential election results by state",
       caption = "Note: Nebraska and Maine split electoral college votes by congressional district") + 
  theme(plot.title=element_text(face="bold", size=12, hjust=.5), 
      plot.caption = element_text(size=10), 
      legend.text = element_text(size = 8))

# Labels
ggplot(stv20, aes(fill = called)) + 
  geom_sf(color = "gray90", linewidth=.2) + 
  scale_fill_manual(values = c("blue", "red"), labels=c("Biden", "Trump")) + 
  theme_void() + 
  coord_sf(expand=FALSE) +
  geom_text(
    aes(label = STUSPS, geometry = geometry),
    stat = "sf_coordinates", size=2, color="white") +
  labs(fill = "",
       title = " 2020 US presidential election results by state",
       caption = "Note: Nebraska and Maine split electoral college votes by congressional district") + 
  theme(plot.title=element_text(face="bold", size=12, hjust=.5), 
        plot.caption = element_text(size=10), 
        legend.text = element_text(size = 8))

# Using "theme_map" (from ggthemes), center title, make a little bigger (expand=FALSE)
ggplot(stv20, aes(fill = called)) + 
  geom_sf(color = "gray90", linewidth=.2) + 
  scale_fill_manual(values = c("blue", "red"), labels=c("Biden", "Trump")) + 
  theme_map() + coord_sf(expand=FALSE) +
  labs(fill = "",
       title = " 2020 US presidential election results by state",
       caption = "Note: Nebraska and Maine split electoral college votes by congressional district") + 
  theme(plot.title=element_text(face="bold", size=12, hjust=.5), 
        plot.caption = element_text(size=10), 
        legend.text = element_text(size = 8))

# Now let's produce a national map at the **counties** level

# First, tell tigris to give us a geographic object at the county level
co <- counties(cb=TRUE, resolution="20m", year="2021", progress_bar=FALSE) |>
  filter(STATEFP != 72) |> 
  shift_geometry()

# How many counties are there in the U.S.?

# Quick view: 
ggplot(data=co) +
  geom_sf(linewidth=.2, fill="white") +
  theme_void()

# We can also overlay this plot with state line borders. Just add another geom_sf command
ggplot(data=co) +
  geom_sf(linewidth=.2, fill="white") +
  geom_sf(data=st, linewidth=.4, color="black", fill=NA) +
  theme_void()

# Read in county-level data from Healy's book (from socviz package)
codata <- county_data

# Match id variables across datasets. We'll match "id" in codata with "GEOID" in tigris
# object.
codata <- codata |> mutate(GEOID=id)

# Merge Healy co. data with geographic maps data from Tigris; will match on GEOID
comerge <- left_join(co, codata)

# Choropleth map of household income at county level. 
# Note use of piping preceding ggplot
# Note: label=dollar is from the "scales" package

# Theme map instead, center title, zooom in a little (expand=FALSE)
comerge |>
  ggplot(aes(fill = hh_income)) + 
  geom_sf(color = "gray90", linewidth=.05) + 
  scale_fill_distiller(palette = "Blues", labels=dollar) + 
  theme_void() + 
  coord_sf(expand=FALSE) +
  labs(fill = "Household 
Income",
       title = "Household Income by County",
       caption = "Source: Healy County Data") + 
  theme(plot.title=element_text(face="bold", size=12, hjust=.5), 
        plot.caption = element_text(size=10), 
        legend.title = element_text(size=10),
        legend.text = element_text(size = 8))

# Reverse legend scale; make high values of income darker colors
comerge |>
  ggplot(aes(fill = hh_income)) + 
  geom_sf(color = "gray90", size=.05) + 
  scale_fill_distiller(palette = "Blues", direction = 1, labels=dollar) + 
  theme_void() + coord_sf(expand=FALSE) +
  labs(fill = "Household 
Income",
       title = "Household Income by County",
       caption = "Source: Healy County Data") + 
  theme(plot.title=element_text(face="bold", size=12, hjust=.5), 
        plot.caption = element_text(size=10), 
        legend.text = element_text(size = 8))

# Overlay state lines on county map; change coloring to yellow-green-blue palette
# Experiment with different palettes and different border colors, etc.
comerge |>
ggplot() + 
  geom_sf(data=comerge, aes(fill = hh_income), color = "gray70", size=.05) + 
  geom_sf(data=st, fill=NA, color="black", linewidth=.1) +
  theme_void(base_size=16) + coord_sf(expand=FALSE) +
  scale_fill_distiller(palette = "YlGnBu", direction = 1, labels=dollar) + 
  labs(fill = "Household 
Income",
       title = "Household Income by County",
       caption = "Source: Healy County Data") + 
  theme(plot.title=element_text(face="bold", size=12, hjust=.5), 
        plot.caption = element_text(size=10), 
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 8))

# Experiment with different palettes and different border colors, etc. Diverging colors
comerge |>
  ggplot() + 
  geom_sf(data=comerge, aes(fill = hh_income), color = "gray70", size=.05) + 
  geom_sf(data=st, fill=NA, color="black", linewidth=.1) +
  theme_void(base_size=16) + coord_sf(expand=FALSE) +
  scale_fill_distiller(palette = "RdYlBu", direction = 1, labels=dollar) + 
  labs(fill = "Household 
Income",
       title = "Household Income by County",
       caption = "Source: Healy County Data") + 
  theme(plot.title=element_text(face="bold", size=12, hjust=.5), 
        plot.caption = element_text(size=10), 
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 8))

# Experiment with different palettes and different border colors, etc.: 
# "scale_gradient2" to define midpoint
comerge |>
  ggplot() + 
  geom_sf(data=comerge, aes(fill = hh_income), color = "gray70", size=.05) + 
  geom_sf(data=st, fill=NA, color="black", linewidth=.1) +
  theme_void(base_size=16) + coord_sf(expand=FALSE) +
  scale_fill_gradient2(low = "red3", mid = "white", high = "dodgerblue2", midpoint = 60000,
                       labels=dollar) + 
  labs(fill = "Household 
Income",
       title = "Household Income by County",
       caption = "Source: Healy County Data") + 
  theme(plot.title=element_text(face="bold", size=12, hjust=.5), 
        plot.caption = element_text(size=10), 
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 8))

#########################################################################################
########################### 2.  World Maps Using Rnaturalearth ##########################
#########################################################################################

# Bring in world map from rnaturalearth
w <- ne_countries(scale = "medium", returnclass = "sf")

# Check "class" of "w." 
class(w)

## [1] "sf"         "data.frame"

#Note "sf" = simple features

# Sort by type and country name
w <- arrange(w, type, name)

# View

# Basic map using defaults; choropleth map for "income_grp" (categorical) variable
# in our "w" object.

w |> 
  ggplot() + 
  geom_sf(aes(fill=income_grp)) + 
  scale_fill_brewer(palette = "YlGnBu")

# Change direction of scale - higher values of income darker
w |> 
  ggplot() + 
  geom_sf(aes(fill=income_grp)) + 
  scale_fill_brewer(palette = "YlGnBu", direction=-1)

# Remove Antarctica and others; remove gray background, relabel legend, theme map
w |> 
  filter(type != "Dependency", type !="Disputed", type != "Indeterminate") |>
  ggplot() + geom_sf(aes(fill=income_grp), color="black", linewidth=.1) + 
  scale_fill_brewer(palette = "YlGnBu", direction=-1, 
                    labels=c("High OECD", "High non-OECD", "Upper middle", "Lower middle", "Low")) +
  coord_sf(expand=FALSE) +
  labs(fill = "Income Level", title = "World Map, Country Income Levels",
       caption="Source: R Natural Earth Data") +
  theme_map() +
  theme(plot.title=element_text(face="bold", size=12, hjust=.5), 
        plot.caption = element_text(size=10), 
        legend.title = element_text(size=10),
        legend.text = element_text(size = 8))

# You can also zoom in using latitudes and longitudes. 

# Projections! Default is crs = 4326; you can change projections

# See:  https://semba-blog.netlify.app/01/26/2020/world-map-and-map-projections/

# For coordinates for specific projections, see: https://epsg.io/

# We're going to use Robinson projection for world maps in this cass: ESRI:54030.
# see NY Times...
w |> 
  filter(type != "Dependency", type !="Disputed", type != "Indeterminate") |>
  ggplot() + geom_sf(aes(fill=income_grp), color="black", linewidth=.1) + 
  scale_fill_brewer(palette = "YlGnBu", direction=-1, 
                    labels=c("High OECD", "High non-OECD", "Upper middle", "Lower middle", "Low")) +
  coord_sf(crs='ESRI:54030', expand=FALSE) +
  labs(fill = "Income Level", title = "World Map, Country Income Levels",
       caption="Source: R Natural Earth Data") +
  theme_void() +
  theme(plot.title=element_text(face="bold", size=12, hjust=.5), 
        plot.caption = element_text(size=10), 
        legend.title = element_text(size=10),
        legend.text = element_text(size = 8))

# Theme map
w |> 
  filter(type != "Dependency", type !="Disputed", type != "Indeterminate") |>
  ggplot() + geom_sf(aes(fill=income_grp), color="black", linewidth=.1) + 
  scale_fill_brewer(palette = "YlGnBu", direction=-1, 
                    labels=c("High OECD", "High non-OECD", "Upper middle", "Lower middle", "Low")) +
  coord_sf(crs='ESRI:54030', expand=FALSE) +
  labs(fill = "Income Level", title = "World Map, Country Income Levels",
       caption="Source: R Natural Earth Data") +
  theme_map() +
  theme(plot.title=element_text(face="bold", size=12, hjust=.5), 
        plot.caption = element_text(size=10), 
        legend.title = element_text(size=10),
        legend.text = element_text(size = 8))

# Let's do a choropleth map for a continuous variable, like GDP
w |> 
  filter(type != "Dependency", type !="Disputed", type != "Indeterminate") |>
  ggplot() + geom_sf(aes(fill=gdp_md_est), color="black", linewidth=.1) + 
  scale_fill_distiller(palette = "YlGnBu", direction=1, labels=dollar) +  
  coord_sf(crs='ESRI:54030', expand=FALSE) +
  labs(fill = "GDP", title = "World Map, GDP",
       caption="Source: R Natural Earth Data") +
  theme_map() +
  theme(plot.title=element_text(face="bold", size=12, hjust=.5), 
        plot.caption = element_text(size=10), 
        legend.title = element_text(size=10),
        legend.text = element_text(size = 8))

# Let's merge outside dataset: VDem (Varieties of Democracy). 
# Let's look at degree of political, rhetorical attacks (by politicians and elites) 
# on the judiciary. 
# Merge in VDem data; I have already created a country id variale ("geounit") that is 
# common between data objects.
v <- read_dta("vdem.dta")

# Merge vdem with world data
vw <- left_join(w, v)

# World map, Robinson
vw |>
  filter(type != "Dependency", type !="Disputed", type != "Indeterminate") |>
  ggplot() + 
  geom_sf(color="black", linewidth=0.1, aes(fill = as.factor(v2jupoatck_ord))) + 
  coord_sf(crs='ESRI:54030', expand=FALSE) +
  scale_fill_brewer(type = "seq", palette = "YlGnBu", na.translate=FALSE,
                    labels=c("Daily/weekly", "Every month", "More than once", "Rare", "None")) + 
  labs(fill = "Attacks", title = "Government Attacks on the Judiciary, 2021",
       subtitle="Based on VDem's v2jupoatck_ord Variable",
       caption="Source: VDem") +
  theme_map() + 
  theme(plot.title=element_text(face="bold", size=15, hjust=.5), 
        plot.subtitle=element_text(size=12, hjust=.5), 
        plot.caption = element_text(size=10), 
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10) )

# Switch direction of scale - darker colors = more attacks
vw |>
  filter(type != "Dependency", type !="Disputed", type != "Indeterminate") |>
  ggplot() + 
  geom_sf(color="black", linewidth=0.1, aes(fill = as.factor(v2jupoatck_ord))) + 
  coord_sf(crs='ESRI:54030', expand=FALSE) +
  scale_fill_brewer(type = "seq", palette = "YlGnBu", na.translate=FALSE,
                    labels=c("Daily/weekly", "Every month", "More than once", "Rare", "None"), 
                    direction = -1) + 
  labs(fill = "Attacks", title = "Government Attacks on the Judiciary, 2021",
       subtitle="Based on VDem's v2jupoatck_ord Variable",
       caption="Source: VDem") +
  theme_map() + 
  theme(plot.title=element_text(face="bold", size=15, hjust=.5), 
        plot.subtitle=element_text(size=12, hjust=.5), 
        plot.caption = element_text(size=10), 
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10) )

# Save as .png file
ggsave("attacks_ord_merc.png", bg="white", w=9, h=4)


###################### ###################### ###################### ###################
###################### Other map options in Tigris (US Maps) #############################
###################### ###################### ###################### ###################

# Zoom in on one state: NY
nyco <- counties("NY", cb = TRUE, progress_bar=FALSE)

# Quick map
nyco |> 
  ggplot() + 
  geom_sf() + 
  theme_void()

# Merge county data to bring in income
nymerge <- left_join(nyco, codata)

# Choropleth map for income for NY counties
nymerge |>
  ggplot() + 
  geom_sf(aes(fill = hh_income), color = "black", linewidth=.1) + 
  theme_void() + 
  scale_fill_distiller(palette = "YlGnBu", direction = 1, labels=dollar) + 
  labs(fill = "Household 
Income",
       title = "Household Income by County, New York State",
       caption = "Source: Healy County Data") + 
  theme(plot.title=element_text(face="bold", size=12, hjust=.5), 
        plot.caption = element_text(size=10), 
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 8))

# Tons of other options in Tigris: counties, cities, census tracts and blocks, 
# school districts
dct <- tracts("DC", cb=TRUE, progress_bar=FALSE)

dct |>
  ggplot() + 
  geom_sf(fill="lightblue", linewidth=.2) + 
  theme_void()

##################################################################################
############################# Other options in RNaturalEarth ######################
##################################################################################

# We can also zoom in on continents; tell rnaturalearth to pull up South America
sa <- ne_countries(scale = "medium", returnclass = "sf", continent = "South America")

# South America
sa |> 
  ggplot() + geom_sf(aes(fill=income_grp), color="black", linewidth=.1) + 
  scale_fill_brewer(palette = "YlGnBu", direction=-1, 
                    labels=c("High OECD", "High non-OECD", "Upper middle", "Lower middle", "Low")) +
  coord_sf(expand=FALSE) +
  labs(fill = "Income Level", title = "Country Income Levels in South America",
       caption="Source: R Natural Earth Data") +
  theme_map() +
  theme(plot.title=element_text(face="bold", size=12, hjust=.5), 
        plot.caption = element_text(size=10), 
        legend.title = element_text(size=10),
        legend.text = element_text(size = 8))

# SA GDP
sa |> 
  ggplot() + geom_sf(aes(fill=gdp_md_est), color="black", linewidth=.1) + 
  scale_fill_distiller(palette = "YlGnBu", direction=1, labels=dollar) + 
  coord_sf(expand=FALSE) +
  labs(fill = "GDP", title = "GDP, South America",
       caption="Source: R Natural Earth Data") +
  theme_map() +
  theme(plot.title=element_text(face="bold", size=12, hjust=.5), 
        plot.caption = element_text(size=10), 
        legend.title = element_text(size=10),
        legend.text = element_text(size = 8))

# South America VDem
vw |>
  filter(continent=="South America") |>
  ggplot() + 
  geom_sf(color="black", linewidth=0.1, aes(fill = as.factor(v2jupoatck_ord))) + 
  coord_sf(expand=FALSE) +
  scale_fill_brewer(type = "seq", palette = "YlGnBu", na.translate=FALSE,
                    labels=c("Daily/weekly", "Every month", "More than once", "Rare", "None"),
                    direction=-1) + 
  labs(fill = "Attacks", title = "Government Attacks on the Judiciary, South America, 2021",
       subtitle="Based on VDem's v2jupoatck_ord Variable",
       caption="Source: VDem") +
  theme_map() + 
  theme(plot.title=element_text(face="bold", size=12, hjust=.5), 
        plot.subtitle=element_text(size=12, hjust=.5), 
        plot.caption = element_text(size=10), 
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10) )

# South America; add country labels
vw |>
  filter(continent=="South America") |>
  ggplot() + 
  geom_sf(color="black", linewidth=0.1, aes(fill = as.factor(v2jupoatck_ord))) + 
  coord_sf(expand=FALSE) +
  scale_fill_brewer(type = "seq", palette = "YlGnBu", na.translate=FALSE,
                    labels=c("Daily/weekly", "Every month", "More than once", "Rare", "None"),
                    direction=-1) + 
  geom_label_repel(
    aes(label = abbrev, geometry = geometry),
    stat = "sf_coordinates",
    min.segment.length = 0,
    size=3,
    point.padding = NA,
    segment.color = "grey20") +
  labs(fill = "Attacks", title = "Government Attacks on the Judiciary, South America, 2021",
       subtitle="Based on VDem's v2jupoatck_ord Variable",
       caption="Source: VDem") +
  theme_map() + 
  theme(plot.title=element_text(face="bold", size=12, hjust=.5), 
        plot.subtitle=element_text(size=12, hjust=.5), 
        plot.caption = element_text(size=10), 
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10) )

## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

#########################################################################################
################### 3.  National and Subnational Maps Using Tigris  ######################
###################    and Tidycensus for Analyzing Census Data    ######################
#########################################################################################

# Go-to source for tidycensus: https://walker-data.com/census-r/index.html
# Analyzing US Census Data: Methods, Maps, and Models in R, by Kyle Walker
# See in particular Chs. 2-6

# To use tidycensus, you need to register for a "key" online from the Census Bureau
# Free and easy to do. The first time you load the package, it will give you instructions
# for getting a key.

# Let's start with educational attainment from 2020 ACS (5 year)
# Go to data.census.gov to look up Census variables. 

# Call up data (and geographical data) from tidycensus 

# National map, subdivided by state; remember, shift geometry for national map, 
# get rid of PR
stateba <- get_acs(
  geography = "state",
  variables = "B06009_005",
  summary_var = "B06009_001",
  year = 2020,
  geometry = TRUE,
  resolution = "20m"
) |>
  shift_geometry() |>
  mutate(propba=estimate/summary_est) |>  
  mutate(pctba=100*estimate/summary_est) |>
  filter(NAME != "Puerto Rico")

# Note: _001 is "total", while _005 is "B.A. degree"

stateba |>
  ggplot(aes(fill = propba)) + 
  geom_sf(color = "black", linewidth=.05) + 
  scale_fill_distiller(palette = "YlGnBu", direction = 1, label=percent) + 
  theme_void() + coord_sf(expand=FALSE) +
  labs(fill = "% B.A.",
       title = "Percent with B.A.",
       caption = "Source: 2020 ACS, 5-year") + 
  theme(plot.title=element_text(face="bold", size=12, hjust = 0.5), 
        plot.caption = element_text(size=10),
        legend.text = element_text(size = 8))

# Interactive map with mapview (see also tmap and leaflet)
mapview(stateba, zcol="pctba")

# County level
coba <- get_acs(
  geography = "county",
  variables = "B06009_005",
  summary_var = "B06009_001",
  year = 2020,
  geometry = TRUE, progress_bar=FALSE,
  resolution = "20m"
) |>
  shift_geometry() |>
  mutate(propba=estimate/summary_est) |>  
  mutate(pctba=100*estimate/summary_est) |>
  filter(GEOID<72000)

coba |>
  ggplot(aes(fill = propba)) + 
  geom_sf(color = "black", linewidth=.05) + 
  geom_sf(data=st, fill=NA, linewidth=.3) +
  scale_fill_distiller(palette = "YlGnBu", direction = 1, label=percent) + 
  theme_void() + coord_sf(expand=FALSE) +
  labs(fill = "% B.A.",
       title = "Percent with B.A.",
       caption = "Source: 2020 ACS, 5-year") + 
  theme(plot.title=element_text(face="bold", size=12, hjust = 0.5), 
        plot.caption = element_text(size=10),
        legend.text = element_text(size = 8))

# Interactive map with mapview (see also tmap and leaflet)
mapview(coba, zcol="pctba")

# Experiment with different palettes. scale_gradient2, define midpoint at average %BA
coba |>
  ggplot(aes(fill = propba)) + 
  geom_sf(color = "black", linewidth=.05) + 
  geom_sf(data=st, fill=NA, linewidth=.3) +
  scale_fill_gradient2(low = "red3", mid = "white", high = "dodgerblue2", midpoint = .20,
                       label=percent) + 
  theme_void() + coord_sf(expand=FALSE) +
  labs(fill = "% B.A.",
       title = "Percent with B.A.",
       caption = "Source: 2020 ACS, 5-year") + 
  theme(plot.title=element_text(face="bold", size=12, hjust = 0.5), 
        plot.caption = element_text(size=10),
        legend.text = element_text(size = 8))

# "School district" level! Let's do MN
mnba <- get_acs(
  geography = "school district (unified)",
  variables = "B06009_005",
  state="MN",
  summary_var = "B06009_001",
  year = 2020, progress_bar=FALSE,
  geometry = TRUE) |>
    mutate(propba=estimate/summary_est) |>
    mutate(pctba=100*estimate/summary_est)
  
mnba |>
  ggplot(aes(fill = propba)) + 
  geom_sf(color = "black", linewidth=.1) + 
  scale_fill_distiller(palette = "YlGnBu", direction = 1, label=percent) + 
  theme_void() + 
  labs(fill = "% B.A.",
       title = "Percent with B.A., Minnesota School Districts",
       caption = "Source: 2020 ACS, 5-year") + 
  theme(plot.title=element_text(face="bold", size=12, hjust=.5), 
        plot.caption = element_text(size=10), 
        legend.text = element_text(size = 8),
        legend.title=element_text(size=10))

mapview(mnba, zcol="pctba")